Importing all the libraries required¶
In [24]:
import fastlmm
import h5py
import numpy as np
import time
# set degree of verbosity (adapt to INFO for more verbose output)
import logging
logging.basicConfig(level=logging.WARNING)
# set figure sizes
import pylab
pylab.rcParams['figure.figsize'] = (10.0, 8.0)
# set display width for pandas data frames
import pandas as pd
pd.set_option('display.width', 1000)
# import the algorithm
from fastlmm.association import single_snp
from pysnptools.snpreader import Pheno, Bed, SnpData
import fastlmm.util.util as flutil
from pysnptools.snpreader import SnpHdf5
from fastlmm.association import single_snp_linreg
#QQ-Plots
from fastlmm.util.stats import plotp
In [3]:
bed_fn = "imputed_snps_binary.hdf5"
pheno_fn = "phenotype.txt"
#Loading the SNP data
f = h5py.File("1001_SNP_MATRIX/imputed_snps_binary.hdf5","r")
#Retrieving the data from the hdf5 snp file.
snps = f["snps"]
accessions = f["accessions"]
positions = f["positions"]
In [4]:
# Details of chromosomes
chr_regions = positions.attrs["chr_regions"]
indices_for_chr1 = chr_regions[0]
print(indices_for_chr1)
In [5]:
#Accessions which are available in 1001genomes org with their respective ids
df = pd.read_csv("AvailableMappedAccessions.csv")
#indices list of accessions that will be used for GWAS
accessions = accessions[:]
accs = accessions.astype(np.int64)
indices_of_accessions = [ np.where(accs == df["Accession Id"][i])[0][0] for i in range(len(df))]
indices_of_accessions = sorted(indices_of_accessions)
print (indices_of_accessions, len(indices_of_accessions))
In [6]:
#Snp data for creating SNP data matrix
#We need iid, sid, and vals
#iid is required to create SNP Data
iid = [[i+1, sorted(df["Accession Name"])[i]] for i in range(len(df))]
print ("iids : ", iid)
genes = snps[:, indices_of_accessions]
#Filtering genes for chromosome 1
genes = genes[: indices_for_chr1[-1]]
threshold= 5
genes = np.transpose(genes)
alleles = np.where(genes.sum(axis=0) >threshold)
snps = genes[:,alleles[0]]
print ("snps shape", snps.shape)
In [7]:
number_of_genes_discarded = genes.shape[1] - snps.shape[1]
print ("number of genes discarded", number_of_genes_discarded)
In [13]:
#sid is required to create SNP data
sid = range(snps.shape[1])
print(sid[:5])
In [14]:
#it's important create snps positions
pos = []
for i in sid:
pos.append([1,i, i])
pos = np.array(pos) #Converting into a numpy array
In [11]:
snps = snps.astype(np.float64)
snps.dtype
Out[11]:
In [15]:
#Creating SNP Data
snp_data = SnpData(iid = iid, sid = sid,
val = snps, pos = pos)
In [16]:
snp_data.sid_count
Out[16]:
In [17]:
snp_data.read().val[0]
Out[17]:
In [18]:
#Verifying Phenotypes
# pheno_vals = pheno_fn.read().val
# print(pheno_vals)
In [26]:
#Working with Phenotypes for features 1-64
phenotype_file_name = "PhenotypesAutoencoder/phenotypesFeature"
suffix = ".txt"
def gwasAnalysisWithManhattanPlot(pheno_fn):
t0 = time.time()
results_dataframe_feature_ = single_snp(test_snps=snp_data, pheno=pheno_fn, count_A1=False)
t1 = time.time()
print ("Time taken for analysis =", str(t1-t0),"seconds")
print("\nManhattan Plot")
flutil.manhattan_plot(results_dataframe_feature_[["Chr", "ChrPos", "PValue"]],pvalue_line=1e-5,xaxis_unit_bp=False)
pylab.show()
print("\nQQ Plot")
plotp.qqplot(results_dataframe_feature_["PValue"].values, xlim=[0,5], ylim=[0,5])
pylab.show()
for feature in range(1,65):
pheno_fn = phenotype_file_name + str(feature) + suffix
# print(pheno_fn)
print("Single SNP GWAS on" +phenotype_file_name+str(feature))
gwasAnalysisWithManhattanPlot(pheno_fn)
